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Abstract 

The software tool GRworkbench is an ongoing project in visual, 
numerical General Relativity at The Australian National University. 
Recently, the numerical diflterential geometric engine of GRworkbench 
has been rewritten using functional programming techniques. By 
allowing functions to be directly represented as program variables 
in C++ code, the functional framework enables the mathematical 
formalism of Differential Geometry to be more closely reflected in 
GRworkbench. The powerful technique of 'automatic differentiation' 
has replaced numerical differentiation of the metric components, re- 
sulting in more accurate derivatives and an order-of-magnitude per- 
formance increase for operations relying on differentiation. 



1 Introduction 

The goal of the ongoing GRworkbench project at The Australian National 
University is to create a visual software tool for numerical operations on 
analytically defined space-times in General Relativity. Such a tool would 
be used by researchers and educators alike to quickly gain insight into the 
physical properties of known exact solutions of the Einstein field equation. 

A new version of GRworkbench was implemented in 1999 [T]. It featured 
a novel numerical differential geometric engine and a flexible visualisation 
system [2 [3], and was easy to extend with additional space-time definitions. 
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Set 


C++ type 


Notes 


TL 
M 
M" 


int 
double 

nvector<double> 

function<B (A)> 


max. ±(2^^ - 1) 

max. ^ ilO'^"^, precision 15 

(as for double) 

see Section [2?T] 



Table 1: Correspondence between certain mathematical sets and C++ 
types in GRworkbench. The term 'max.' refers to the largest representable 
elements of the set; the term 'precision' refers to the number of significant 
figures to which elements of the set are represented. 



In 2003, the numerical and differential geometric aspects of GRwork- 
bench were rewritten using functional programming techniques, enabling 
the direct representation in the CH — h code of GRworkbench of those con- 
cepts in numerical computation and Differential Geometry that are nor- 
mally defined in terms of functions. The new functional framework en- 
ables potentially complex numerical experiments to be quickly and sim- 
ply defined [1]. Numerical experiments performed in GRworkbench were 
employed in our analysis of a recent scientific claim, described elsewhere 

Functional programming and its applications in GRworkbench are de- 
scribed in ij2l Automatic differentiation, in which the derivatives of a func- 
tion are exactly computed whenever the value of the function is computed, 
is described in ij3l In f|4] we outline planned future developments for the 
numerical engine of GRworkbench, in which interval arithmetic will be used 
to establish a guaranteed bound on numerical errors. 

2 Functional programming 

In the traditional programming languages of scientific computing, programs 
typically consist of routines that operate on data stored in program vari- 
ables. Every variable in C++ has a type, and there is an approximate 
correspondence between C++ types and standard mathematical sets. Ta- 
ble [T] lists some important examples. 

The first two sets in Tableware, of course, directly representable in some 
way in every language of scientific computing. The nvector<T> type uses 
the C++ template mechanism ([5], page 327) to provide a type representing 
n-tuples of any other type T. Thus, nvector<double> represents elements 
of M" and nvector<nvector<double>> represents elements of M.mxn, the 
set of m X n matrices with real- valued entries. 

^For the original scientific claim see [?]■ 
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The following is a routine in C++: 



double mean(double a, double b) 

{ 



return (a + b) / 2; 

} 



The corresponding mathematical definition is 



mean: K x M 



R, 



mean(a, b) 



a + b 



(1) 
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The first line of the routine conveys the same information as the first line 
of HI): the routine mean takes two real numbers as arguments, and returns 
a real number. The first line of the routine without the routine name 
or argument names, i.e., double (double, double), is the signature of the 
routine; it conveys the same information as R x K ^ R in ([T|). The rest of 
the routine definition, enclosed in braces, encodes the second line of ([T]). 

We may define a function as anything which behaves like the routine 
mean above, in the sense that it accepts zero or more arguments, and 
returns a value. In the most common programming languages of scientific 
computation, C and Fortran, the only possible functions are routines, and 
so the terms 'function' and 'routine' are used interchangeably. The key 
feature of functional programming is that there can be functions other 
than the routines typed in by the programmer — functions created while the 
program is running. The mechanism to achieve this in C++ is introduced 
in tJ2.1.11 we first describe how functions can be stored in variables in C++. 

2.1 Functions as data 

The capability to store functions in variables is not unique to functional 
programming. Most languages used for scientific computation have some 
way to store a reference to a program routine. For example, in C and 
C++ the address &if of a routine double f (double x) can be stored in a 
function pointer variable of type double (*)(double). Function pointers 
can be called just like routines. 

GRworkbench uses the Boost Function library [9] to store references 
to functions. The Boost Function library provides the templatised type 
function<T> representing a function whose signature is T. The following 
code fragment shows how the routine mean can thus be stored in a variableo 

function<double (double, double)> f = mean; 
// the following two lines are now equivalent 

^Anything after the characters // in a Hne of code is ignored by the CH — h compiler. 
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double X = f(l, 2); 
double X — mean(l, 2); 

Observe from the last two lines that the variable f can be used just like the 
routine mean; they are both functions. 

In general, if we let (^i x • • • x An B) denote the set of functions from 
Ai X ■ ■ ■ X An to iJ, then the corresponding C++ type in GRworkhench 
is function<B (Ai, . . . , An)>, where the sets B, Ai, . . . , A„ correspond to 
the types B, Ai, . . . , Ap. The fourth row of Table [T] summarises this rela- 
tionship. 

The most important application of the storage of functions in variables 
is that functions can then be arguments to other functions. Consider the 
following routine, which crudely approximates the derivative of a function 
/ at a point x: 

double slope(function<double (double)> f, double x) 

{ 

double h = 0.1; 

return (f(x + h) - f(x - h)) / (2 * h); 

} 

The corresponding mathematical definition is 
slope: (M ^ M) X R ^ M, 

, N f{x + h)- f{x-h) , „ (2) 
s\o^eU,x)^^^ '-^ h = Q.l. 

Many other numerical algorithms naturally take a function as an argu- 
ment. Two examples are 

minimise: (K ^ R) x R ^ M, 

(3) 

minimise(/, x) — (a local mininram of / near a;), 

and 

integrate: (M ^ M) x R x R ^ M, 



integrate(/, a, h) — (numerical estimate of / f{x) dx). 

J a 



(4) 



2.1.1 Creating functions at run-time 

Consider the following function, defined in terms of the slope function ([2]): 
derivative: (R ^ R) ^ (R ^ R), 

(5) 

derivative(/) = (7:R— i-M, g{x) — slope{f,x). 
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Concept 


Representation in GRworkhench 


Coordinates 
Metric components 
Inter-chart map 
Point 

Tangent vector 

Metric 

World-line 


nvector<double> 
nvector<nvector<double>> 
see Section |2.2.2| 
point 

tangent_vector 

function<double (tangent_vector, tangent_vector)> 

function< point (double)> 



Table 2: Representation of important differential geometric concepts in 
GRworkhench. 



For each function /, it returns the function which returns the slope of f at 
its argument. 

In most traditional languages of scientific computing it is not possible 
to encode a routine which returns derivative(/) for all functions / : R ^ K. 
In C-I--I- it is possible to encode this function, by means of a functor class 
([8], pages 514-515). The result is a C++ routine function<double (double 
)> derivative(function<double (double)> f). (For complete code, including 
that of the necessary functor class, see [B], pages 13-14, 92-94.) 

If we were to replace the slope routine with a more sophisticated al- 
gorithm for numerical differentiation, then this derivative routine would 
be a good approximation to the mathematical operation of differentiation. 
For example, derivative(sin) would be a good approximation to the function 
cosH Functional programming permits numerical operations, like derivative, 
to be expressed in a way which closely resembles the mathematical opera- 
tions that they approximate. 

2.2 Functional Differential Geometry 

Many fundamental notions in Differential Geometry and General Relativ- 
ity, such as the action of the metric tensor, and particle world-lines, are 
functions. By elevating functions to the same level as traditional data types 
(Z, R) , functional programming makes these notions directly representable 
as variables in the C++ code of GRworkhench. Table [2] summarises the 
correspondence between important concepts in Differential Geometry and 
their representations in GRworkhench. 

•^Many standard functions, including sin and cos, are built-in to C++. 
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2.2.1 Charts and the metric components 

On the n-dimensional space-time manifold A4, a chart is a pair (U,4>) 
representing a coordinate system on the set U C M, where the one-to-one 
function (f>: U ^ M" maps points in U to their coordinates in R". Space- 
times are defined in GRworkbench by the specification of the components 
of the metric tensor on (f>{U) C K." for one or more charts (U,(j>), and 
by the definition of the maps 0^ o between the coordinate systems of 
overlapping pairs of those charts, where o denotes function composition. 

The coordinates of a point on a chart, a;* € K", are represented by 
a variable of type nvector<double>. The components gab of the metric 
tensor at a point on a chart are represented as an n x n matrix, by a 
variable of type nvector<nvector<double>>. A function which defines the 
metric components as a function of the chart coordinates might then be 
of the form 

chart: ^ M„xn, 

(6) 

chart(x') = 5afcU', 

represented in GRworkbench by a function of signature nvector<nvector< 
double>> (nvector<double>). In general, however, the chart coordinates 
are an open subset of M", and so ^ will not be defined everywhere in M". 
A mechanism is required to represent functions which are only defined on a 
subset of some other, standard set. (By 'standard set' we mean a set which 
is already represented by a type in C-I--I-, such as those listed in Tables [T] 
andU) 

GRworkbench employs the Boost Optional library [lOj to represent func- 
tions which are undefined for some values of their arguments. The Boost 
Optional library provides a templatised type optional<T>, which repre- 
sents the set S U {0}, where S is the set corresponding to the template 
parameter type T, and is a special value taken by functions at points 
where they are undefined. 

Using the optional mechanism we rewrite ([6|) to support charts defined 
on subsets of R": 

chart: R" M„xri U {0}, 

/ i\ I 9ab\x^j if the are valid chart coordinates; (7) 
chart (a; ) = < 

I 0, otherwise. 

The corresponding C++ type is 

function <optional<nvector<nvector<double>>> (nvector<double>)>, (8) 

for which GRworkbench declares the short synonym chart, using the C-l — h 
typedef mechanism: 
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typedef function<optional<nvector<nvector<double>>> (nvector< 
double>)> chart; 

The optional mechanism is most useful when the caller of a function 
cannot know beforehand whether the function will be defined at the ar- 
guments to be given to it. This is the case for callers of functions of the 
form ([7]). The optional mechanism thus enables the differential geometric 
and visualisation algorithms in GRworkbench to be coded in such a way 
that they can operate on any space-time definition, without prior knowl- 
edge of the particular coordinate systems (charts) in which they will be 
working. 

2.2.2 Inter-chart maps 

For two charts (Ua,(t>a) and {Lip, (f)f3), the inter-chart map from <f>a{Ua) to 
(t>(3{Ui3) is 

where </'/3|w„ is the function (pp restricted to the set Ua- 

The domain of (pap is, in general, a subset of R". Hence (pap cannot 
be represented by a variable of type function<nvector<double> (nvector< 
double>)>; instead, the optional mechanism is again employed. Thus, an 
inter-chart map from a chart {lAa, pa) to a chart {Up, 4>p) is represented by 
a function 

map: M" ^ K" U {0}, 

'('/'/3k if [x') e MUc) and cp-\x')€Up; 

0, otherwise. 

(10) 

The corresponding C-I-+ type is 

function<optional<nvector<double>> (nvector<double>)>. (11) 

The C++ typedef mechanism is used to define a synonym map for this 
type. 



map (a;*) 



2.2.3 World-lines 

The abstract notion of a point p S A^, independent of any particular 
coordinate system, is represented in GRworkbench by a C-l — I- type point. A 
point is constructed from two pieces of information: a chart which contains 
it, and its coordinates on that chart. 
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The operator[ ] routine of the point type takes one argument, a variable 
of type chart, and returns a variable of type optional<nvector<double>>, 

representing the coordinates of the point on the given chart. (The optional 
mechanism is used because a particular point may not have coordinates on 
the given chart.) Thus, if p is a variable of type point, and c is a variable of 
type chart, then the coordinates of p on c are given by p[c]. If there is no 
inter-chart map defined from the point's original chart to the chart c then 
is returned. 

A general curve in space-time, such as a world-line, which may not 
be defined for all values of its parameter, is a function A: K — »■ U {0}; 
such functions are represented by variables of type function<optional<point 
> (double) >. Curves in GRworkbench must be represented in this form, 
because they are often defined in terms of numerical algorithms that may 
not converge everywhere to a solution, even if a solution exists; such a 
curve returns the value for parameter values for which there was no 
convergence to a solution. The synonym worldline is defined for the type 
function<optional<point> (double)>. 



2.2.4 Tangent vectors 

The abstract notion of a tangent vector v G Tp, where Tp is the tangent 
space of a point p E A4, is represented in GRworkbench by the C-l — h 
type tangent_vector. A tangent_vector is constructed from three pieces of 
information: the point to whose tangent space it belongs, a chart containing 
that point, and the contravariant components of the tangent vector on that 
chart. 

As with the point type, the operator[ ] routine of the tangent_vector 

type takes one argument, a variable of type chart, and returns the com- 
ponents of the tangent vector on the given chart, in a variable of type 
optional<nvector<double>>. When the components of a tangent vector 
are requested on a chart other than that from which the tangent vector was 
constructed, GRworkbench uses the inter-chart map, if it exists, to com- 
pute the components. If are the components of a tangent vector w at a 
point p on a chart with coordinates x*, then the components on another 
chart, with coordinates , are 



- v\ (12) 



The entries of the matrix A\ are the derivatives of the inter-chart map 
(p: M" M" with respect to the coordinates a;* of its argument, evaluated 
at p. In its most recent version, GRworkbench computes A\ , and thereby 
the components , using the method of [J3l 
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At a point p, the metric gab is naturally considered as the inner product 

metric : T„ x r„ ^ R, 

\ (13) 
metric(u, w) — gabu'^v . 

This function is encoded in GRworkhench in the routine metric, whose 
signature is double (tangent_vector, tangent_vector). Also, the operator* 

routine of the tangent_vector type is defined to call metric, so that if u and 
V are variables of type tangent_vector, then the expression u * v is equivalent 
to the expression metric(u, v). This notation is reminiscent of the two 
equivalent forms 

gabu'^v'' = Ubv'' (14) 
for the inner product of two tangent vectors. 

2.3 Numerical operations and applications 

We may think of the operation of numerically determining a geodesic by 
integrating the geodesic equation from initial conditions as a function 

geodesic: Tp ^ (K ^ X U {0}), 
geodesic(w) = A, A:M^AlU{0}, 

where A is the numerically determined geodesic with tangent vector v at 
p = A(0). Similarly, the numerical parallel transport of a tangent vector 
along a curve may be represented by a function 

paralleLtransport : (R ^ 7W U {0}) x T{M) ^ (R ^ T{M) U {0}), 
paralleLtransport(/, v) ^ h, /i : R T{M) U {0}, 
h{t) is the parallel transport of w G 7/(o) to f{t) along /, 

(16) 

or h{t) = if, for example, the numerical parallel transport algorithm did 
not converge to a solution. 

In terms of these two functional definitions it is trivial to define the 
following interesting object: 

paralleLcurve: (R ^ X U {0}) x T{M) ^ (R ^ X U {0}), 
paralleLcurve(/, w) =7, 7: R U {0}, (17) 

7(t) = geodesic(paralleLtransport(/, z;)(i))(l), 

which represents the world-line of an object which is stationary at a proper 
distance ygabV^v^ with respect to an observer whose world- line is /. By 
'gluing' together functional objects such as these, it is easy to define poten- 
tially complex numerical experiments simulating interesting physical situ- 
ations [5]. 
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3 Automatic differentiation 



Differentiation plays a role in most numerical operations in GRworkbench, 
because the derivatives of the components of the metric tensor feature 
in, for example, the geodesic equation and the parallel transport equation. 
Previously in GRworkbench, numerical differentiation of arbitrary functions 
f : W V, where V is any vector space, was accomplished via numerical 
estimation of the limit 

lim f(- + h)-fi--h) ^ 

using the technique of Richardson extrapolation ([5], pages 18-21; [11| . 
pages 186-189). In (jlSp . d{h) is the centred difference approximation to 
the derivative of / at x. 

While this method is easily applicable to any function, it has some im- 
portant drawbacks. Firstly, it requires many evaluations of the function /, 
at X ± h for numerous values of h. Secondly, an estimate of the limit (|18p 
may not converge, or, what is worse, may converge to a value that is incor- 
rect. The accurate convergence of the algorithm depends on the values of 
h at which d{h) is evaluated being of approximately the same size as the 
smallest scale over which the function varies significantly around x. For a 
general implementation of the numerical differentiation method (|18[) . this 
information may not be available. Thirdly, the accuracy of this method is 
rarely more precise than half as many significant figures as the precision 
of the floating point arithmetic employed. (For the double type this is 
~ 15/2 ~ 7 significant figures.) 

The technique of automatic differentiation ([12j. Chapter 5) avoids all of 
these difficulties. The technique follows from three observations. The first 
is that all C-f-l- routines encoded by the programmer must, at the lowest 
level, eventually be defined in terms of a finite set of built-in fundamental 
operations (such as addition, extraction of square roots, and trigonometric 
functions). The second is that all of these fundamental operations have ex- 
act derivatives at (almost) every point at which they are defined, and these 
exact derivatives may themselves be expressed in terms of the fundamental 
operations. The third is that from the exact derivatives of the fundamen- 
tal operations we can obtain the exact derivative of any function defined 
in terms of the fundamental operations, through the use of the chain rule 
for differentiation: 



^f{u{x)) = ^f{u) 
dx du 



^u{x). (19) 

u{x) 



These facts are exploited in the method of forward automatic differ- 
entiation, which is implemented in GRworkbench. In this method, the 
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fundamental CH — h operations which act on real numbers are extended to 
also act on ordered pairs {u, u') of real numbers, with u being the usual ar- 
gument to the function, and u' representing the derivative of u with respect 
to some independent variable x. For example, we redefine 



In this way, the derivative of a function is automatically computed 
whenever the value of the function is computed. The derivative is accurate 
up to round-off error in the floating point arithmetic, and its accuracy is 
independent of the behaviour of the function around x. As well as being 
more accurate, the method is also faster than numerical differentiation, 
because the extra computation involved in computing the second part of 
each ordered pair {u,u') (roughly a factor of two) is generally less compu- 
tationally expensive than the many evaluations of / that are necessary to 
numerically estimate the value of the limit . 

The replacement of numerical differentiation by automatic differentia- 
tion has resulted in a large performance and accuracy increase for those 
numerical operations in GRworkbench that are dependent on the differen- 
tiation of the components of the metric tensor. 

4 Future directions 

Currently in GRworkbench, tangent vectors and metrics, respectively ten- 
sors of type (1,0) and type (0,2), are each implemented separately. We 
will replace them with an implementation of general tensors of arbitrary 
type, whose components on any chart are determined in terms of their 
components on one chart by the general transformation rules. This im- 
plementation will be based around a C-I--I- type that can represent any 
function of vectors and covectors that is linear in each of its arguments. 
Thus, the representation of tensors in GRworkbench will be nearly identical 
to their usual mathematical definition. 

4.1 Interval arithmetic 

Many computations require us to make statements not about points, but 
rather about open sets of points. Ordinary differential equation solvers 



sin: R X M ^ R X M, 
sm{u,u') — (sin M, u' cos It) 



(20) 



and 



multipUcation: (M x R) x (R x R) ^ M x M, 
(m, u'){v, v') = {uv, u'v + uv'). 



(21) 
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require that functions be defined, bounded, and have some number of con- 
tinuous derivatives not at a single point, but in a neighbourhood around 
an initial condition. Just as we can implement automatic differentiation by 
extending the definition of basic operations, we can also implement interval 
arithmetic (|12j. Chapter 3), so that for 

/:M^M (22) 

we define a new function 

/':MxM^KxM (23) 

that satisfies 

/'(a, 6) ^{c,d)^yxe [a,b]J{x) G [c,d]. (24) 

The function /' is a function from one interval of real numbers to another 
interval of real numbers. By combining the definition of a function in 
terms of fundamental operations with a bound on its parameter, we obtain 
a bound on the value of the function. For example, the function 

multiply : M X M ^ M, 

, , (25) 
multiply(a, c) — ac, 

is generalised to act on intervals of real numbers by 

multiply' ; (M X M) X (R X M) ^ (M X R), 

multiply' ((a, b), (c, d)) = (min(ac, ad, be, bd), max(ac, ad, be, bd)) . 

Thus, the product of the interval (2.9, 3.1) with the interval (3.9, 4.1) is the 
interval (11.31,12.71). 

When combined with automatic differentiation and the optional mech- 
anism, a whole range of statements about any function defined in terms of 
fundamental operations can be tested on any interval. We might be able 
to establish with certainty that a function is everywhere defined on a cer- 
tain closed 'box' [a,b] x [c,d] x .. . C R", that it is bounded, or even that 
it is C''. These properties are important to guarantee that geodesic trac- 
ing and other algorithms can respond correctly to singularities and other 
pathological conditions. 

A concrete example occurs for the metric components of the Schwarz- 
schild space-time expressed in standard spherical polar coordinates, which 
are valid separately for < r < 2M and for 2M < r < oo. Presently in 
GRworkbench, two charts are used, respectively covering r € (0, 2M) and 
r € (2M, oo), so that a numerical operation, such as geodesic tracing, that 
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docs not 'know' about the r = 2M singularity, cannot accidentally 'skip 
over' it by discretely sampling points on either side of it. 

Using interval arithmetic it would, in principle, be possible to employ 
one 'pseudo-chart' covering the whole range < r < oc, because any nu- 
merical operation that evaluated the metric components on an interval 
containing r = 2M would learn that the metric components are unde- 
fined somewhere in that interval, and are unbounded over it the algorithm 
should then attempt to investigate a smaller subset of the original interval. 

A more widely applicable further benefit of interval arithmetic is that 
it can provide powerful control over round-off errors. The exact result for 
any calculation typically lies between two representable floating-point num- 
bers; round-off error occurs when one of these must be chosen to represent 
the result, a process analogous to writing a value to only so many decimal 
places. Instead of choosing one value, an interval arithmetic operation re- 
turns both values, which together bound the tightest representable interval 
enclosing the true result. Interval arithmetic trades the illusory precision 
of floating point computations for guaranteed accuracy. 

Once it is implemented in GRworkbench, the 'pseudo-algebraic' method 
of interval arithmetic, which exploits exactly-known properties of the fun- 
damental mathematical operations in C-|— 1-, will enable a powerful new 
method of exploring properties of space-times and render the ability to 
make exact statements about these properties. 

5 Conclusion 

The new functional framework for the differential geometric engine of GR- 
workbench represents the underlying mathematical structure of Differen- 
tial Geometry more closely than ever before. Functional programming 
makes it easier to define systems that model interesting physical situa- 
tions in numerical experiments. The technique of automatic differentiation 
provides accurate and fast derivatives for functions, such as analytically- 
defined space-time metrics, that are defined in terms of certain fundamental 
mathematical operations. In the future, interval arithmetic may enable the 
numerical engine of GRworkbench to be written in such a way that there 
are no unknown numerical errors at all. 
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